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Abstract 

The first three dynamic multipole polarizabilities for the ground state of hydrogen, helium, 
hydride ion, and positronium hydride PsH have been computed using the variational Monte 
Carlo (VMC) method and explicitly correlated wave functions. Results for the static dipole 
polarizability by means of the diffusion Monte Carlo (DMC) method and the finite field ap- 
proach show the VMC results to be quite accurate. From these dynamic polarizabilities van der 
Waals dispersion coefficients for the interaction of PsH with ordinary electronic systems can be 
computed, allowing one to predict the dispersion energy for the interaction between PsH and 
less exotic atoms and molecules. 

PACS number(s): 36.10.-k, 02.70.Lq 
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While experimentalists relay everyday on positrons and positronium atoms (Ps) to collect informa- 
tion about microscopic features of macroscopic systems like solutions, polymers and crystals, much 
less effort has been devoted to the theoretical understanding of the complex interactions that take 
place between ordinary matter and positrons. Among the explored avenues of this field, we mention 
the interest in predicting the stability of classes of compounds like e + M and MPs [1-9], where M 
represents an atomic or molecular system, and the calculation of the cross sections in the scattering 
process of e + and Ps on a molecule or an atom [10-16]. 

On the contrary, the evaluation of the interaction energy between e + M or MPs and a molecule 
or atom is an almost unexplored issue fTijj . We believe this fact is primarily due to the need of 
a very accurate trial wave function to describe correctly the correlated motions of electrons and 
positrons. So far, only variational calculations with explicitly correlated Gaussians (ECG) jj], || or 
Hylleraas-type functions [7-9], and the DMC method [2-4] have shown to be able to adequately 
recover the correlation energy in positron-containing systems. 

Related to the calculation of the interaction energies is the calculation of second order properties 
of positron-containing systems, a problem whose surface has been barely scratched in the past 
H. These properties, specifically the dynamic polarizabilities, are strictly related to the van der 



Waals coefficients that describe the long range interaction between systems 18 , representing a way 
to tackle the problem of the asymptotic intermolecular interactions. Recently, Caffarel and Hess 
showed that these properties can be computed by means of quantum Monte Carlo simulations []l9| 
connecting the imaginary-time-dependent dynamics of the unperturbed system with the transition 
probabilities of a reference diffusion process. In this work we apply a modified version of their method 
to compute dynamic multipole polarizabilities for PsH, H, He, and H~ as a way to understand the 
behaviour of these systems when interacting with an external field, and as a first step towards the 
definition of the interaction potential between PsH and the ordinary matter. 

As far as we know, the work by Le Sech and Silvi jTJ] is the only one reporting calculations on 
the effect of a constant electric field on PsH. In that work they computed both the static dipole 
polarizability, 123 a.u., and the behaviour of the annihilation rate T2 7 versus the intensity of the field 
employing explicitly correlated wave functions, numerical integration, and a variation-perturbation 
approach. As by-product of our calculations of the potential energy curve of the e + LiH system j2TJ , 
we obtained an estimation of the static dipole polarizability of 49(2) a.u., a value quite different from 
the one computed by Le Sech and Silvi. Since we believe this difference to be too large to admit an 
explanation based on the different accuracy of the methods used to compute this value, we plan to 
solve this puzzle in this work. 

In the method by Caffarel and Hess [|l9) the frequency-dependent second order correction to the 
ground state energy is written as a sum of the two time-centered autocorrelation functions of the 
perturbing potential V 

/>oo 

E { 1\oj) = - e ±tul C vv {t)dt (1) 
Jo 

where the autocorrelation function Cyv (t) is given by 

Cw(t) = (V(0)V(t)) 9 , - {V(0)\ n (2) 

Here, indicates that the average has to be taken using the Langevin dynamics that samples 

the square of the exact ground state wave function of the unperturbed system. 

Caffarel and Hess |l!| showed that it is possible to compute Cvv(t) employing an optimized 
trial wave function and the pure-diffusion Monte Carlo (PDMC) method, an alternative algorithm 
to the commonly used DMC with branching, where each walker explicitly carries its own weight 
along all the simulation [^l| . 

In their work on He and H2, Caffarel, Rerat, and Pouchan (22) reported that the autocorrelation 
function Cvv(t) becomes dominated by the noise at large times, and this fact might be due to the 
fluctuations of the walker weights that increase during a PDMC simulation, while the value of the 
autocorrelation function itself becomes smaller. While the second effect is intrinsic to the stochastic 
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method, the first can be reduced employing a more accurate trial wave function that is able to 
reduce the weight fluctuations. Another possibility, giving up the exactness of the method (i.e. not 
sampling the exact \&o), is represented by the sampling of a quite accurate trial wave function without 
carrying around the weight for each walker, a method we call Perturbation Theory variational Monte 
Carlo (PT-VMC). This algorithm can be useful for those systems whose autocorrelation function 
has a large decaying time, as the case of H~ and PsH. This large decaying time will increase the 
fluctuations of the carried weights, and hence the statistical noise in the autocorrelation functions 
in the long-time region. 

As a test of the correctness of our computer program and of the accuracy of the method, we 
computed the first three autocorrelation functions, and hence the dynamic polarizabilitics up to the 
octupolar one, for the two systems H and He. The analytical forms of the perturbing potentials were 
taken from Ref . |22f . While for H we employed the exact ground state wave function and compared 



with the analytic values of the multipole polarizability 18 , for the He case we used a 25 term 



Hylleraas-type wave function optimized by means of the standard energy minimization [|23| . We 
fitted the numerical Cvv(t) results of our simulations with a linear combination of three exponential 
functions 
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C vv (t) ~J2<He~ Xit (3) 

i=l 

in order to have an analytical representation of the autocorrelation functions at all the times. Since 
it is important to reproduce accurately the long time behaviour of Cvv, the smallest \i in Eq. [| 
was independently calculated fitting In [Cvv] i n the long time region with a first order polynomial. 
This choice was found to improve sensibly the goodness of the total fitting in this time range. 

These analytical representations of Cvv allow us to compute easily the integrals in Eq. [I] and to 
obtain simple expressions of a{oj). The parameters obtained by the fitting procedure are available 
from the authors upon requests. 

For both systems we found excellent agreement of the static polarizabilities (H adip =4.495 au, 
OLquad =15.034 au, a oc t =133.105 au; He adip =1.382 au, a qua d =2.401 au, a oc t =10.367 au) with 
the exact results for H p8[ , with PDMC results by Caffarel, Rerat, and Pouchan p2| , with Glover 
and Weinhold upper and lower bounds for He |24|] , and with the accurate results by Yan et al. |25| . 

At this point we would like to stress that, although in the PT-VMC method the walkers carry 
always a unitary weight because the branching process is absent, similarly to the PDMC method 
the time step has to be chosen short enough to produce only a small time step bias. For these 
two systems we found the time step of 0.01 hartree -1 to be adequate to compute statistically exact 
results. 

As a check of the ability of the PT-VMC method to compute polarizabilities also for highly 
polarizable systems whose exact wave function is more diffuse than the one of He and H, we selected 
the hydride ion as test case. For this system we optimized a 5 term Hylleraas-type wave function 
whose average properties are shown in Table |l] together with the accurate results obtained in Ref. 
]26[ ] . Table [l] contains also the multipole static polarizabilities computed in this work employing a 
time step of 0.01 hartree -1 , and the static polarizability computed by Glover and Weinhold ]24| . 
Comparing the mean values in Table 0, one can notice that our 5 term wave function gives lower 
values than the ones obtained in Ref. ]26[ | except for (r_ ) . This fact may explain the underestimation 
of the adip by PT-VMC, that recovers 92(2) percent of the accurate value. Nevertheless, this result 
represents a fairly good estimation of the static dipole polarizability for H , a quantity that appears 
difficult to compute even with more complex approaches |27l . 

As far as PsH is concerned, we computed the autocorrelation functions using two different trial 
wave functions, including one (\&y.) and 28 (\&§?) terms ||. The choice of two trial wave functions 
to guide the Langevin dynamics was aimed at testing the dependency of Cvv(t) on the quality of 
the wave function itself. 

Employing the PT-VMC method and our wave functions for PsH, we computed the autocorre- 
lation functions for three perturbation potentials: 
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Vi = x\ + x 2 - x. 



(4) 
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where the subscripts 1 and 2 indicate the two electrons, while the subscript p indicates the positron. 
These potentials are the cartesian forms of the dipole, quadrupole, and octupole moment operators 
for the PsH system. Figure 1 shows the averaged correlation functions for Vi, Vi, and V3 as obtained 
by the VMC method employing the 28 term trial wave function. Each value of the correlation 
functions was computed employing roughly 10 10 configurations. From Figure 1 one can note the effect 
at large evolution times of the dispersion of the "trajectories" used to compute the autocorrelation 
function. This effect makes difficult the reproduction of the long-time regime of these functions due to 
the exponential decay and the roughly constant statistical error introduced by the method. Moreover, 
the statistical error strongly depends on the perturbation potential whose autocorrelation function 
is computed, i.e., more specifically on the dispersion of its mean value over the ^fj, distribution. 

The results for the static multipole polarizabilities, i.e. for to — 0, computed with both trial 
wave functions, are shown in Table ||. While for the dipole polarizabilities there is a good agreement 
between the two values, larger differences are present for the higher multipole polarizabilities. This 
fact is an indication of the different accuracy of the two functions in approximating the exact wave 
function at large distances from the nucleus. In fact it can be shown that if one approximates 
the autocorrelation functions taking care only of the excitation to the first state of the appropriate 
symmetry, the autocorrelation function is proportional to (V 2 ) — (Vi) 2 , where Vi is the perturbing 
potential. Comparing the dipole results with the value obtained by Le Sech and Silvi ||, again 
the large difference of the computed polarizabilities is apparent. As a final check for this problem, 
we computed the energy of the PsH when immersed in a weak static electric field F by means of 
standard DMC simulations adding the linear potential F[x\ + x 2 — x p ). To make our simulations 
stable, i.e. to avoid the dissociation of the PsH, we truncated the effect of the linear potential at 

= 15 bohr. We fitted the DMC results by means of the simple polynomial a — a.di p F 2 /2, where 
ctdip is the static dipole polarizability, obtaining aa v — 42.3(8) a.u. We believe that this result, 
statistically indistinguishable from the a dip obtained by the PT-VMC method, gives the definitive 
answer to the problem of the PsH polarizability. Nevertheless, the discrepancy between our PT- 
VMC and DMC adi P and the one computed by Le Sech and Silvi |j) remains puzzling. In our 
experience to compute the matrix elements they needed, millions of configurations must be 
used even for systems like PsH to avoid to be fooled by a false convergence. Unfortunately, Le Sech 
and Silvi did not report any information about the number of configurations they used to compute 
the integrals, so we cannot judge the numerical accuracy of their results. 

An attempt to estimate the total accuracy of our a results can be made comparing the polariz- 
ability values obtained by the two wave functions. These differ by 10 percent at most, a value that 
we feel might give a conservative estimate of the relative errors for the higher multipolar fields. 

As stated previously, although dynamical polarizabilities are interesting on their own, they repre- 
sent the basis to compute van der Waals dispersion coefficients for the interaction between different 
systems. Therefore, following Ref. (25), we present the calculation of the Cq, Cg, and C10 dispersion 
coefficients between H, He, and PsH as a first effort to obtain accurate information on the interac- 
tion between positronic systems and ordinary matter in the framework of the Born-Oppcnheimer 
approximation and second-order perturbation theory. 

Using the fitted parameters for H, He, and PsH we computed the coefficients for the interaction 
between the ordinary systems and between these and PsH. The values are reported in Table ^. Since 
the values for the H-H, H-He, and He-He coefficients are accurately known pq ], we use them as a 
test of the accuracy of our approach: all the values differ from the accurate results by Yan et al. 
|25[ ] at most by one part over hundreds. 
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Comparing the C n 's for the ordinary systems with the ones for the interaction with PsH, it strikes 
that these last are more than an order of magnitude larger than the formers. These features, due 
to the larger PsH polarizability, indicates that positronic systems strongly interact with ordinary 
matter even at large distances. Unfortunately, nothing can be said about location and depth of 
the total potential minimum. This strongly depends also on the effect of the repulsion between the 
positron cloud and the H and He nuclei, so that we believe a supermolecule approach is needed. 
In a previous work |l7| we computed the interaction energy between H and PsH, showing that 
this system could have a metastable state. Although the dispersion coefficients for the interaction 
between He and PsH are smaller than those for PsH and H, they might be large enough to give rise 
to a potential well that could support at least a stable state. If this turns out to be the case, the 
He-PsH system could be the lightest van der Waals (i.e. bound by means of dispersion forces) stable 
dimer. 
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Figure captions: 

Figure 1. Logarithm of the correlation functions of the perturbing potentials. 
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Table 1: Mean values for observables of the ground state 1 5' of H . All values are in atomic units. 



This work (5 term wave function) 



b Ref. 
c Ref. 
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876.9(3) 


34848(71) 



Table 2: Static multipole polarizabilities for the ground state 2,1 S of the PsH computed with one 
term an d 28 term (^y 8 ) wave functions. All values are in atomic units. 
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Table 3: Computed dispersion coefficients. All values are in atomic units. 

Ref. H 
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